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Method and Apparatus for quantifying tissue histology 

This inyention relates to a method and apparatus for quantifying tissue 
histology. In particular the invention relates to methods using an analysis of 
the spectra of remitted light to establish information on the properties of the 
tissue. The invention is also applicable to methods and apparatus which rely 
upon a spectral analysis of light remitted, emitted, and/or transmitted from 
any material or object under test where they have parameters exhibiting 
wavelength specific optical effects. 

There exists the need for a system, which can recover histological 
parameters from biological tissue in a way which is invariant to the 
intensity of the incident illumination and scene geometry. It is an objective 
of the present invention to provide such a technique. Such a system would 
be of value in systems where the topology of the tissue or image surface in 
not known a priori. It would also be of value in. a system where the intensity 
of the illuminating light cannot be assumed constant. Potential applications 
include but are not limited to imaging and analysis of the tissue of the 
gastrointestinal track with an endoscope and imaging and analysis of skin 
over areas where there is a significant change in curvature, such as the face. 

A system is currently in existence, which is able to assist clinicians in their 
diagnosis of melanoma. The technique is based on a patent, international 
patent application publication number W098/22023. This system is based 
25 on the discovery that when the range of colouration of normal human skin 
is plotted in a standard RGB colour space, it lies on a well-defined surface. 
Furthermore, if an abnormality such as dermal melanin is present, the 
colouration of the skin changes in such a way as to move points away from 
the surface which describes healthy skin. By incorporating a calibration 
30 calculation which allows variation of dermal thickness to be taken into 
account, the technique is able to detect abnormalities and thus assist 
clinicians in their diagnosis of melanoma. 
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The fundamental principle behind this system is that it is possible to 
construct a mathematical function that relates image values, measured using 
a digital camera, to appropriate histological parameters. Using this 
5 . functional relation, it is possible to obtain the value of each parameter at 
every point across a given image. A parametric map can then be produced 
which gives a grey-scale representation of the parameter value across the 
whole image. 

10 Although this system has been proved to be clinically effective, it requires 
exact calibration of the illuminating light source and does not take into 
account any variation in surface geometry. Thus the technique is limited to 
problems where a probe can be placed in contact with the region of interest. 
This ensures that the incident light is controlled and calibrated and that it*s 

15 angle of incidence remains constant. 

The proposed invention relates to a method for imaging tissue in such a way 
to give quantitative spectral data independently of the surface geometry or 
the intensity of the illuminating light. This will allow a non-contact form of 
20 imaging and analysis which will be applicable to many different 
applications. The method may be used with the technique described in 
W098/22023 and subsequent related patents but is not exclusive to it. 

The method concentrates upon the analysis of light remitted by the tissue, ie 
25 the illuminating light which penetrates the tissue to some depth and is 
reflected (or scattered or/and absorbed) to different degrees at different 
depths due to different parameters of the tissue. Effects due to surface 
reflection are to be eliminated from the analysis, 

30 Substantial work has been carried out to develop image analysis algorithms 
which are able to identify different objects irrespective of the illuminating 
light. Many of the techniques developed are based around the linear model 
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of surface reflectance as proposed in L. Maloney and B. Wandell, "Color 
constancy: a method for recovering surface spectral reflectance", J. Opt. 
Soc. Am. A 3, 29-33 (1986). This approach is based on the idea that the 
surface reflectance of any object within an imaged scene can be expressed 
5 as a weighted sum of basis spectral reflectance functions: 

SiX) = ±ajSj{X) (1) 
i=i 

and that the illuminating light can similarly be expressed as a weighted sum 
of basis lights. It has been shown that only a small number of basis 
functions are required to obtain accurate approximations to the surface 
10 reflectances of many naturally occurring objects and also the spectral 
variation of natural daylight. 



With this technique it is possible to recover the vector of weighting 
constants aj from a vector of image values and thus specify the spectral 

15 reflectance of the imaged object at every pixel. Every potential imaged 
object/object characteristic will have a unique spectral reflectance. Thus, if 
the spectral reflectance can be determined using a linear model, then the 
parameter vector can be specified. With this approach it should be possible 
to recover a parameter vector from the vector of image values at each pixel. 

20 Unfortunately the method is only able to recover the weighting constants aj 
to within a multiplicative scaling factor and thus cannot be used to specify 
the exact spectral reflectance and therefore the exact parameter vector. 



An approach to geometry-insensitive segmentation of images has been 
25 developed in G. Healey, "Using colour for geometry-insensitive 
segmentation," J. Opt. Soc. Am. A 6, 920-937 (1989), and is based on the 
idea of normalised colour. With this approach image values are first 
divided by an estimate of normalised colour. This estimate is based on 
approximating the incoming signal from colour pixel values by using a 
30 finite-dimensional linear approximation to represent the colour signal. 
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Using these normalised values, different metal and dielectric materials can 
be identified across an imaged scene in which the geometry varies 
considerably. 

5 A similar technique has been applied to evaluate bum injuries {M. A. 
Afromowitz, G. S. van Liew and D. M. Heimbach, "Clinical evaluation of 
bum injuries using an optical reflectance technique", IEEE trans. Biomed. 
Eng. BME-34, 114-127 (1987), and M. A. Afromowitz, J. B. Callis D. M, 
Heimbach, L. A. Desoto and M. K. Norton, "Mulitspectral imaging of bum 
10 wounds: a new clinical instrument for evaluating bum depth", IEEE tran. 
Biomed. Eng. 35, 842-849 (1988)}. In this case, RGB image values were 
normalised by dividing them by the response of an IR filter. From the 
normalised values it was possible to assess the extent of bum damage 
across a given area of imaged skin. 

15 

There exists a need for a non-invasive technique for analysing an object or 
material (which may be complex, for example multi-component and/or 
multi-layer and which may be solid, gaseous, liquid, etc) which does not 
require calibration to take into account changing illumination conditions, 

20 

Summarv of the Inventions 

According to a first aspect of the invention there is provided a method of 
analyzing at least one parameter of a body component, comprising the steps 

25 of illuminating the component with light of at least a first and second 
waveband, receiving light of at least said first and second wavebands 
remitted by the component at a photoreceptor, and analyzing the light 
received at the photoreceptor to provide a ratio between the amount of light 
remitted of the first waveband and the amount of light remitted of the 

30 second waveband , and from this calculating the component parameter. 
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Thus the invention lies in the appreciation that by skilful selection of the 
wavebands of light remitted by a biological component, usually human or 
animal tissue, the ratio between two such wavebands can be used to create a 
useful parametric image of the biological component . The wavebands may 

5 be calculated using a biological or mathematical model to monitor the 
relationship between a particular image ratio and the parameter to create a 
function which then can be used for monitoring that parameter in the 
biological component. As an alternative to creating a function the measured 
waveband ratios can be compared with the predictions of a model either 

10 mathematical or experimentially measured. 

According to a second aspect of the invention, there is provided a method 
of analyzing at least one parameter of a body component, comprising the 
steps of illuminating the component with light of at least a first and second 
predetennined waveband, receiving light of at least said first and second 
predetermined wavebands remitted by the component reflected by the 
surface at a photoreceptor but eliminating light reflected by the surface of 
the component, where the predetermined wavebands are chosen such that 
the component parameter is a one to one function of the ratio between the 
amount of light remitted by the component of the first predetennined 
wavebands and the amount of light remitted by the component of the second 
predetermined waveband, and analyzing the light received at the 
photoreceptor to provide a ratio between the light of the first waveband and 
the light of the second waveband , and from this calculating the component 
parameter. 

The method is such that the effects of reflection rather than remittance will 
be ignored. Although this would appear to limit the application to 
components which do not have a specular component of reflection, such as 
30 many organic objects, image processing algorithms have been developed to 
allow removal of this component of reflection giving greater applicability 
of the technique. For example, it has been shown that it is possible to 
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remove the highlights from complex images containing inhomogeneous 
dielectrics. It is also possible to remove the surface component of 
reflection using a polarising filter. Once this component has been removed 
from image data, it will be possible to use the techniques described here. 



Preferably the effects of surface reflection are eliminated by providing a 
pair of cross-polarised linear polarizing filters. The first of these is placed 
in front of the source of illumination and the second in front of the image 
capture system. There are however other methods which will be apparent 
10 the skilled reader which could be used to elinunate surface reflection 
effects. 

The body component may be any biological component but is most usefully 
animal tissue. 



Each waveband referred to may comprise a single wavelength, but in 
practice will preferably comprise a band of wavelengths, detectable by the 
photoreceptor. 

20 The light emitted by the light source may be a white light or light of a 
plurality of wavelengths, some of which are outside the predetermined 
wavebands, and filters may be used to limit the light received by the 
photoreceptor to the desired wavebands. Such filters may be placed between 
the light source and the tissue, between the tissue and the photoreceptor or 

25 at both positions. Alternatively white light may be emitted by the light 
source and received by the photoreceptor with the analysis means 
establishing the amount of light within the desired waveband. 

To understand why this process eliminates any variation in illumination 
30 intensity and surface geometry it is necessary to consider the dichromatic 
reflection model. This was first proposed by S. Shafer in "Using colour to 
separate reflection components" Color. Res. Appl. 4, 210-218 (1985) and 
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States that light remitted from an object is the sum of two components, the 
"^body" component and the surface" component. The body component 
refers to physical processes occurring after penetration of light into the 
material and the surface term to reflections which take place at the surface 
5 of the object. The body component is a function of the spectral 
characteristics of the object, whereas the surface component depends only 
on the object geometry and the incident light. The theory states further that 
each component can be considered the product of a geometrical term and a 
wavelength dependent term. 

10 

The proposed invention is used where an optical system makes spectral 
measures of tissue. One embodiment of such an optical system uses a 
colour digital camera as the photoreceptor, although a monochrome digital 
camera arranged to take sequential images with different coloured light 

15 sources could also be used. Both these embodiments may have millions of 
image pixels or very few pixels, or even just one pixel in the case of the 
monochrome system. The optical system may work in the visual spectrum, 
or over an extended spectrum to include non visible wavelengths. These 
non visible wavelengths may include infra-red light. This infra-red light 

20 may include wavelengths in the 600nm to SOOnm band. 

In the case of a conventional colour digital camera, the system measures 
light through a number of optical filters. Image values for a specific image 
location, corresponding to the nth filter, are given by 

25 

where Kb and Ks are the geometric terms of the body and surface 
component respectively and Cb and Cs are colour terms. By using the system 
30 of polarising filters described above it is possible to eliminate surface 
reflection. Image values are then given as a simple product of a geometric 
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term and a colour, or wavelength dependent term. The illuminating light is 
now written as 

where f^is a wavelength independent scaling factor determined by the 

intensity of the light source but which does not change, or chnages in a 
known manner, will wavelength. This allows the dichromatic reflection 
model to be written as 



i" = € jE^iA)SiA)R\A)dA 



where e — s^K^. The function R^(X) defines the spectral response of the nth 
filter and S^(A) the remitted spectrum of the illuminated tissue. It is 
15 essential that both jBq ('^) R^'C^) are known for the given imaging 
system. Thus the invention is preferentially utilised in to systems where 
tissue of interest is illuminated with light of know spectral characteristics. 

If the optical system records a M-dimensional vector of image values at 
20 each pixel then it is possible to define a iV-dimensional vector of image 
ratios, which is obtained by defining appropriate ratios of image values. An 
example of such a vector is 




As the constants depends only on position within an image all components 
of the ratio vector r will be independent of the constants and thus 
independent of the illumination intensity and geometrical factors in the 
imaged scene. 
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The invention is applicable to problems in which all histological variation 
can be described by K parameters. The concept of a parameter vector is 
introduced and defined as 



where the space P defines all possible parameter variation and thus 
variation in tissue histology. Using the current invention it is possible to 
recover a parameter vector from a vector of image ratios. To achieve this it 
is necessary to have some technique for predicting a vector of image ratios 
from a given parameter vector. This can be achieved via some experimental 
technique or with an appropriate mathematical model of light transport 
within the tissue of interest. Techniques such as Monte Carlo modelling or 
the Kubelka-Munk approximation have been developed for this purpose. 
With such a model it is possible to predict a remittance spectrum which 
corresponds to a unique point in parameter space, that is a unique tissue 
histology. With a knowledge of the spectral response of the illuminating 
light source and the spectral response of the filters, used in the image 
acquisition system, it is possible to predict a vector of image values for a 
given point in parameter space. This can be expressed as 



where the space I defines all possible measurements made by the optical 
system. Using an appropriate definition of image ratios, such as on the one 
given above, it is possible to obtain a vector of image ratios. This can be 
expressed as 



P = {Pn/^2v...Px) 




re R 
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where the space R defines all possible iitnage ratios that can be obtained 
from the space of image measurements. A function / can now be defined 
which maps from points in parameter space to points in the space of image 
ratios. To implement this function it is first necessary to compute the 
5 spectral reflectance of the material of interest for the given set of parameter 
values, or point in parameter space. Using this spectral reflectance, along 
with the spectral responses each of the filters a vector of image 

values can be calculated. Finally from this a vector of ratios can be 
obtained. This three-stage mapping can written as 

10 

to denote the mapping from parameters space to the space of image ratios. 
Provided that a remittance spectrum can be defined for any possible 
15 parameter combination then this mapping is defined for the whole or 
parameter space. The proposed invention deals with the inverse of this 
function, defined as 

20 

which denotes the mapping from the space of image ratios back to 
parameter space. A key part of the invention is to establish whether a 
suitable function g can be defined which will allow any measured ratio to 
be mapped back to the appropriate parameter combination. Such a mapping 

25 must be 1-1. That is, for every point in the space of image ratios there must 
be a corresponding unique point in parameter space. If this is not the case, 
ambiguity will arise as it could be possible to recover more that one set of 
parameter values from a given vector of image ratios. To establish this 
condition, it is first necessary to deal with the function /, which must be 

30 considered a vector valued function of a vector variable, that is. 
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r = /(p). 



To establish whether this function is 1-1 the determinant of the Jacobian 
matrix, corresponding to this mapping, can be analysed. This is defined as 



J = 



3/, a/i 3/i 




3pi dp2 dpfc J 



If the determinant of this matrix is non-zero at a point in parameter space 
then there exists a neighbourhood around this point where the function /can 
10 be approximated linearly. This means that any points within this region will 
map under a 1-1 mapping to a unique point in parameter space. If, when 
using a system to image a given tissue, it can be established that the 
Jacobian is non-zero across the whole of parameter space then the function/ 
will be 1-1 everywhere. 

15 

Once this condition has been established it is necessary to find either an 
approximation or an exact analytic expression for the function g which will 
enable image ratios to be mapped to specific parameters. Although in some 
cases it may be possible to obtain an analytic function, in most cases it will 

20 be necessary to construct a piecewise continuous approximation. This can 
be achieved by discretising parameter space in suitably small intervals and 
generating the corresponding image ratio values for every point within a 
discretised space. Some form of multidimensional interpolation technique, 
such as a cubic spline, is then used to construct a continuous piecewise 

25 approximation to the function g. This then allows processing of pixels from 
an imaged tissue to give the corresponding set of parameter values. Any 
number of pixels may be processed in this way to produce a number of 
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parametric maps, which give quantitative information on the parameters of 
interest across the whole of the imaged scene. 

Such maps are of immense value to clinicians and other persons interested 
5 in the composition of specific tissues. 

The implementation of the proposed invention proceeds along the following 
steps: 

10 1. For the tissue to be imaged identify all parameters whose variation 

could cause a change in spectral remittance when illuminated with 
light. 

2. Have, by some means, a method for predicting the spectral 
remittance of a given tissue for any combination of the identified 

15 tissue parameters. 

3. Establish the spectral responses of each channel of the given imaging 
system and from this define an appropriate set of image ratios. 

4. Check that the mapping from the space of parameters to the space of 
image ratios is 1-1 over the range of all parameter variation. 

20 5. If this condition holds obtain some function, either exact or 

approximate, which maps points in the space of image ratios to the 
corresponding point in parameter space. 
6. Using this function images can then be processed to give quantitative 
information on the underlying tissue histology. 

25 

According to a third aspect of the invention there is provided apparatus for 
analyzing at least one parameter of a body component, comprising a light 
source for illuminating the component with light of at least a first and 
second predetermined waveband, a photoreceptor for receiving light of at 
30 least said first and second predetermined wavebands remitted by the 
component reflected by the surface at a photoreceptor; surface reflection 
elimination means for eliminating light reflected by the surface of the 
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component, where the predetermined wavebands are chosen such that the 
component parameter is a one to one function of the ratio between the 
amount of light remitted by the component of the first predetermined 
wavebands and the amount of light remitted by the component of the second 
5 predetermined waveband, and microprocessor means for analyzing the light 
received at the photoreceptor to provide a ratio between the light of the first 
waveband and the light of the second waveband, and from this calculating 
the component parameter. 

10 The wavebands having image ratios which map 1-1 to a parameter of 
component vary depending upon the particular component, and the 
particular parameter to be analysed. 

Typically the method and apparatus are used to analyse all the parameters 
15 required to characterize a particular component, with the light source and 
photoreceptor emitting and receiving for each parameter, a pair of 
wavebands chosen such that the ratio between the amounts of light remitted 
by the component of each waveband (ie the image ratio for that pair of 
wavebands) is a 1-1 function of the particular parameter. In practice, the 
20 minimum number of wavebands to be monitored will be equal to n+1, 
where n equals the number of parameters. 

It has been found for skin that three parameters characterize the tissue, 
namely skin thickness, melanin concentration and blood concentration and 
25 melanin and blood concentration may be analysed effectively using the 
methods and apparatus of the invention. 



The required predetermined wavebands may be found using the method 
described above iteratively. 
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According to a fourth aspect of the invention, there is provided a method of 
deriving a pair of predetermined wavebands suitable for use in analysing a 
given parameter of a body component, the method comprising the steps of : 
1) defining a set of potential wavebands 
5 2) defining one or more image ratios, the or each image ratio for a 

region being obtained by dividing the amount of light remitted by the 
component of a given waveband for that region, the "image value" 
for that filter, by another image value for that same region; 

3) for a component parameter to be analysed and for said defined set 
10 of potential wavebands and for said given image ratios, obtaining a function 

mapping points in parameter space to points in image ratio space; 

4) determining whether the mapping function provides a 1:1 
correspondence between points in parameter space and points in image ratio 
space; and 

15 5) if the mapping function does not provide a 1:1 correspondence, 

rejecting said potential wavebands, repeating steps 1) to 4) and, if the 
mapping function does provide a 1:1 correspondence accepting the potential 
wavebands as a candidate set of predetermined wavebands. 

,20 The key step in the present invention is that of identifying a set of filter 
properties and image ratio and a mapping function which maps image ratios 
to tissue parameters with a 1:1 correspondence. This first step may require 
consideration of many potential filter properties and corresponding mapping 
functions. 

25 

Preferably, for each filter the method of the present invention is used to 
determine the centre wavelength of the filter. The method may additionally 
be used to determine the full width half maximum (FWHM) of the filter 
characteristic. 

30 

Preferably, step 3) comprises constructing a Jacobian matrix for the 
mapping function with respect to said parameter(s), and obtaining the 
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determinant of that matrix. If the determinant is strictly positive or strictly 
negative over the entire parameter space, then a 1:1 correspondence 
between points in parameter space and points in image ratio space is 
assumed. It will be appreciated that these operations may be performed 
5 using multiple processing steps or may be combined into a single 
processing step. 



Embodiments of the present invention may involve the calculation of an 
error, indicative of the accuracy of parameter recovery obtained using said 
10 mapping function. The error may be calculated as follows: 

a) calculate the error associated with image acquisition for each 
vector 

of each image ratio; 

b) from the image ratio vector error, calculate the maximum 
15 possible error in each component of the parameter vector across the 

whole of parameter space; and 

c) use the vector of parameter errors at each point within parameter 
space to measure the accuracy of parameter recovery. 



20 Alternatively corrections may be made for error by standard mathematical 
error correction algorithms the choice of which will be apparent to the 
skilled addressee of the specification. 

The present invention may be incorporated into many different filter 
property calculation schemes. For example, in a scheme using a genetic 
25 algorithm, the method may be used to identify a plurality of candidate filter 
parameter sets. The method of the present invention is then applied 
repeatedly to find an optimal filter parameter set using the candidates. 

Alternatively, the method may be employed in a scheme using a gradient 
30 descent algorithm. In such a scheme, the method of the third aspect of the 
present invention is employed to identify a first candidate set of filter 
parameters. This set is then used to make a selection of another set of filter 
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properties, and the process repeated as often as necessary to arrive at an 
optimal solution. Of course, any suitable optimization algorithm can be 
used to compute an optimal solution or a solution which has sufficient 
accuracy. 

5 

Each time the method steps 1) to 4) are carried out, the image ratios may or 
may not be changed. That is to say that, for each repetition, the potential 
wavebands and the image ratios may be changed, or only the potential 
wavebands may be changed. 

10 

Although this invention is applicable with particular advantage to the non 
invasive analyzing of tissue - typically animal and preferably human tissue 
it will be appreciated that the method and apparatus could also be used to 
monitor parameters of a material where the parameters characterizing the 
15 material have wavelength specific optical properties and where it is 
possible to control the spectral characteristics of the illumination. 



According to a fifth aspect of the invention there is provided a method of 
20 determining a property or properties of each of a set of filters, which filters 
are used to select specific wavelength ranges to quantify a parameter or 
parameters of a tissue, the method comprising the steps of: 

1) defining a set of potential filter properties 

2) defining one or more image ratios, the or each image ratio for a 
25 region being obtained by dividing the quantified output of a given filter for 

that region, the "image value" for that filter, by another image value for 
that same region; 

3) for an object or material to be analysed and for said defined set of 
potential filter properties and for said given image ratio, obtaining a 

30 function mapping points in parameter space to points in image ratio space; 
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4) determining whether the mapping function provides a 1:1 
correspondence between points in parameter space and points in image ratio 
space; and 

5) if the mapping function does not provide a 1:1 correspondence, 
5 rejecting said potential filter properties, repeating steps 1) to 4) and, if the 

mapping function does provide a 1:1 correspondence accepting the potential 
filter properties as a candidate set of filter properties, 

10 According to a sixth aspect of the present invention there is provided 
apparatus for analysing an object or material having means for conducting a 
spectral analysis of remitted, emitted, and/or transmitted light to quantify a 
parameter or parameters of an object or material, the apparatus comprising 
a plurality of filters for splitting said light into respective components, the 

15 filters having properties obtained by using the method of the above first 
aspect of the invention. 

It will be appreciated that the filters of the apparatus may be implemented 
in optical, electrical/electronic, or software form, 

20 

Brief Description of the Drawings 

Methods and apparatus according to the various aspects of the invention 
will now be described, by way of example only with reference to the 
25 accompanying drawings, in which:- 

Figure 1 illustrates a function f which maps an image ratio i to a material 
parameter p; 

Figure 2 illustrates a combination of errors in parameter space, mapped to 
30 image ratio vector space; 

Figure 3 illustrates a model of the layered structure of skin; 
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Figures 4a and 4b show how the remitted spectrum (intensity vs 
wavelength) varies for different melanin and blood levels respectively; 
Figure 5 illustrates a set of filters suitable for analysing blood and melanin 

levels in skin; 

5 Figure 6 is a flow diagram illustrating a method of defining a set of filter 

properties for use in analysing the properties of an object or material; 

Figure 7 is a flow diagram illustrating a method of defining a set of suitable 

wavebands for use in analyzing the parameters of tissue; 

Figure 8 is a schematic view of a method and apparatus for analyzing facial 
10 skin; and. 

Figures 9a , 9b and 9c are respectively, a colour image of a human face 
taken with a standard digital camera and parametric maps, showing a grey- 
scale representation of then quantitative measurements of melanin and 
blood derived using the method and apparatus in accordance with the 
15 second and third aspect of the invention.. 

The proof of the theory behind the selection of appropriate wavebands and 
image ratios for a given parameter will now be described with reference to 
figures 1 ,2 and 6 

20 

In a typical analysis system, light remitted from an object can be measured 
using a digital camera with a small number of optical filters representing a 
number of wavebands. Image values - brightness or "intensity" - for each 
image location (x, y) for a given filter (the nth filter) are given by: 

= Kt f E{X)S{X)R'{X) dX + Ks f E(X)Br{X) dX (2) 
25 J J 

where and Ks are the geometric terms of the body and surface 
components respectively and and are colour terms. The first integral 
in equation (2) is the product of three terms: E(A) is the illuminating light, 
S(A) is the spectral remittance from the body of the imaged object, and 
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is the spectral response of the nth optical filter. In the second integral 
there are only two terms as there is no wavelength dependence on the 

surface component of reflection. The dichromatic reflection model is very 
important for 3-D scene analysis as it allows for both colour and 
5 geometrical analysis of objects within a scene. 



A key issue is to show that the technique proposed here is valid for 
problems where the intensity of the illuminating light is unknown (whilst 
assuming that the spectral definition of the illuminating light is known). 
10 For this purpose the incident light is written as: 

E{X) = €ojBo(A) (3) 

where eo is a wavelength independent scaling factor. Equation (2) now 
becomes 

i^(x, y) = ej J7o(A)5(A)ii^(A) dX (4) 

15 where e = 8oKb. A digital camera records an iV-dimensional vector of image 
values at each location (x; y). If a mapping, which is independent of the 
constant 8, can be established between the vector of image values and the 
vector of parameters, then it will be possible to recover scene parameters 
from image data in a way that does not depend on illumination intensity or 

20 scene geometry. 

We now introduce the concept of an image ratio, obtained by dividing one 
image value, calculated from equation (4), by another. For a given image 
vector, the nth image ratio is given as: 



25 



^ = ^ n^m. (5) 

Simple consideration of equation (4) shows that any ratio defined in this 
way will be invariant to a change in the parameter e. Thus any method for 
the recovery of parameter values from image ratios will be independent of 
scene geometry and illumination intensity. 
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The objective here is to extract quantitative parameters upon which the 
object colouration depends, not to find statistical similarities. Moreover, 
the specific filters are chosen to maximise the distance in the image ratio 
5 space between vectors corresponding to similar parameter values, as this 
minimises the error on the parameter value recovered from the colour 
image. 

The technique described is generally applicable to scenes in which a small 
10 number of parameters are required to describe all possible objects/object 
characteristics. In the formulation, the parameters will be considered to 
vary continuously. Thus, the technique will be particularly applicable to 
problems where object characteristics need to be measured across an image. 
For example a medical imaging system may be required to analyse a 
15 particular tissue. The underlying structure of the tissue will not vary, only 
specific characteristics such as thickness of the different layers (including 
zero thickness) or the concentration of a particular chemical constituents 
(including zero concentration). In this situation a small parameter vector 
can describe all possible variations in the characteristics of the imaged 
20 scene. For K scene parameters the parameter vector is defined as: 

K 

P = El^* P^I" (6) 

and the space P defines all potential object characteristics. Ultimately, a 
mapping from image ratios back to the parameter vector is required, but 
first the forward problem of obtaining image ratios for a given parameter 
25 vector is considered. A reflectance spectrum, corresponding to a given 
point within parameter space, can be described by the vector in M 
dimensional wavelength space: 

Af 

A = 5^ Am A € A 

m=l ^ ' ^ 

where the space A defines all possible spectral reflectance functions. The 
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mapping a, defined as 

is introduced to denote the mapping from parameter space to wavelength 
space. This mapping gives the spectral reflectance of the object specified 

5 by the vector p. Such a mapping can be achieved either by a spectroscopic 
measurement, or by using a mathematical model which takes as input the 
parameters and produces a corresponding reflectance spectrum. Models of 
light propagation in different media, such as the Monte Carlo method or the 
Kubelka Munk approximation, can be used for this purpose. It must be 

10 possible to perform this mapping across the whole of parameter space, thus 
defining every possible spectral reflectance function. 



15 



20 



A digital camera with N optical filters records an iV-dimensional image 
vector at each pixel. The image vector is given as: 

i=f:i" i€i (9) 

where I describes the space of all possible image values. The process of 
image acquisition can be considered as the projection of points in 
wavelength space to points in filter space. This projection is performed by 
the mapping function: 



b:A 



(10) 



Equation (4) performs this mapping b in continuous form. In discrete form, 
the response of the nth optical filter, is given as: 

M 

^^ = ^E^A„, (11) 

m=l 

where Am= Eo(/l)S(/l) and the positive weights at each wavelength are given 
25 by , thus defining each filter response function. A digital camera 
effectively performs the mapping b, projecting points from a large 
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dimensional space (wavelength space) to a small dimensional space (filter 
space). With such a mapping there will be a substantial loss of information. 
However, even with this loss of information, it should be possible to define 
the mapping in such a way that accurate information regarding the original 
5 parameter values can still be recovered from image data. Conditions for 
this mapping will be discussed in the following later. 

Most current image acquisition systems use an RGB system of filters. 
Although this defines a potential mapping b, it may not be the best mapping 

10 with which to recover parameter values from image data. However, it is 
known to select specific filters to obtain better clarity of data than that 
possible with an RGB system (although mainly for visualization or image 
segmentation, not for quantification). Also, in spectrometry, particular 
spectral wavelengths are selected using statistical methods to improve 

15 quantification of components in mixtures. It will therefore be appreciated 
that an objectively defined set of optical filters is able to perform the task 
of recovery of parameters, which describe the variation in human skin, 
better than a standard RGB system. 



20 Once the vector of image values has been obtained, a vector of image ratios 
can be 

calculated using equation (5). The vector of image ratios is given as: 

ivr 

where I describes the space of all possible image ratios. The mapping from 
25 filter space to the space of image ratios is performed after image acquisition 
and will be referred to as mapping c, defined as: 

c:I— >I. (13) 

There are many ways to define the image ratios and thus the mapping c. For 
example, pairs of image values could be used to define image ratios as: 
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20 



n=l,2,...,— (14 

or a single image value could be taken as the denominator with which to 
calculate image ratios from the remaining image values, for example: 

?^ = ^ n = 2,3,...,JV. (15) 
^1 

5 At most the dimensionality, N , of the new space will be one less that that 
of the original filter space, This would correspond to the definition 
given in equation (15). Alternatively, if the image ratios were defined as 
given in equation (14), then the dimensionality of the new space will be 
half that of the original filter space. The aim is to recover a J?-dimensional 

10 parameter vector from N image ratios. Thus there must be at least as many 
image ratios as parameters, that is, N >K . 

The function /defined as: 

/ = ao&oc /:P — 5^1 (16) 

15 represents the three stage mapping from parameter space to wavelength 
space, to image space, and finally to the space of image ratios. For a given 
set of optical filters, it will be possible to perform this mapping across the 
whole of parameter space, provided that it is possible to obtain a spectrum 
for any given parameter vector. The inverse of function /is defined as: 



/-I : I p (17) 



and maps from the space of image ratios directly to parameter space. If it is 
possible to define an appropriate / '\ it will be possible to recover 
parameter values from image data in a way that is independent of 
illumination intensity and scene geometry. The ultimate aim is to find the 
25 optimum / which maximises the accuracy of parameter recovery. Before 
a detailed discussion of this mapping is presented, it is important to 
emphasise that the form of the function / will depend on the mappings a, b 
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and c. Although mapping a is fixed for a given problem, mapping b will 
vary with the choice of optical filters and mapping c will vary depending on 
how the image ratios are defined. 



5 Any mapping function which is to map from the space of image ratios (I- 
space) to parameter space (P-space) must be 1 to 1. That is, for a given 
point in P-space, there must be a corresponding unique point in I -space and 
vice-versa. If this is not the case, ambiguity will arise as it could be 
possible to recover more that one set of parameter values from a given 

10 vector of image ratios. Once this condition has been established, it is 
necessary to consider the error' associated with parameter recovery as, using 
a digital camera, it will only be possible to obtain image values to within a 
given uncertainty. This will introduce an uncertainty into the final recovery 
of the parameter vector. There could also be an error associated with the 

15 prediction/measurement of a spectrum from the parameter vector. For 
simplicity the analysis presented here will be restricted to problems in 
which the error associated with the spectral measurement can be neglected. 

Initially, the problem where one parameter is sufficient to describe all 
20 variation in an imaged scene will be analysed. The methodology will then 
be extended to problems where the number of parameters is greater than 
one. 



Consider the case where one image ratio (two image values) is used to 
25 recover a single parameter value. Figure 1 illustrates a function / which 
gives the image ratio as a function of the parameter p. It is clear that in 
order to satisfy the 1 to 1 condition, the curve must not have any turning 
points: that is, it must increase or decrease monotonically in the appropriate 
range of p. Mathematically this is expressed as: 



30 



dp 



> 0 Vp € P. 



(18) 
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Measurement of an image ratio value Iq, corresponding to a parameter value 
po, is now considered. Associated with acquisition of each image value is 
an uncertainty due to camera error. It is straightforward to show, using 
standard error analysis, that the error associated with an image ratio i , 
5 which has been calculated from the two image values ij and 12, is given as: 



where Ai is the camera uncertainty. This error has been shown on the 
ordinate of the graph in Figure 1. If the derivate of /is non-zero in some 



10 linearly. Assuming the error Af to lie within this neighbourhood, the 
corresponding error in the parameter value is given as: 



Thus, it is possible to obtain a value for the error Ap, associated with 
parameter recovery, at any point in P-space. An optimisation criterion can 

15 then be defined based on some measure of this error. For most applications 
it will be necessary to minimise the error equally across the whole of P- 
space. For others it may be that high accuracy parameter recovery is 
required within a certain range of parameter values. For example, in a 
medical image application, imaged tissue could be deemed pathological 

20 once a characterising parameter changes beyond a threshold level. This 
would need to be accounted for with some form of weighting in the 
optimisation criterion. It is interesting to note that in order to niinimse Ap, 
it is necessary to maximise the magnitude of the derivative given in 
equation (18). This will ensure that any search, carried out to find an 

25 optimum /, will tend to move towards regions of search space where the 1 
to 1 condition is satisfied. 




(19) 



neighbourhood of po then it is possible to approximate this function 




(20) 



In theory it is possible to recover the parameter using more than one image 
ratio. In this case it will be necessary to calculate the error associated with 
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parameter recovery for each of the image ratios and select the one, at each 
point in P-space, which has the smallest associated error (Ap). It may be 
that the optimisation procedure gives a single image ratio which performs 
better than any other across the whole of P-space. In this situation there is 
5 no benefit to using more that one image ratio. 

The analysis is now extended to the general problem where the recovery of 
a iiT-dimensional parameter vector is required from an N dimensional 
vector of image ratios. Initially the analysis will be restricted to the case 
10 where N =^ K and will then be extended to include situations where N > K. 
As discussed earlier, if iV < jRT, then it is not possible to recover K- 
dimensional data from an N dimensional measurement. 



15 



20 



25 



The mapping function f , defined as: 
i = f(p) 



(21) 



must now be considered a vector valued function of a vector variable. In 
the following analysis specific results from differential geometry will be 
used. For further details the reader is directed to for example M. M. 
Lipschutz, Differential geometry (McGraw-Hill Book Company, New York, 
1969), To establish whether the function /provides a 1 to 1 relationship, it 
is first necessary to consider the behaviour of the determinant of the 
Jacobian matrix, simply referred to as the Jacobian. This is defined as: 



det 



\dpk) 



dp2 



dpk 



(22) 



The Jacobian can be considered the multidimensional equivalent of the one 
dimensional derivative given in equation (18). The inverse function theorem 
states that, if the Jacobian is non-zero at a point po in P-space, then there 



wo 2004/010862 ^■pCT/GB2003/003367 



27 

exists a neighbourhood around po where the function f can be approximated 
linearly as 

f (p) = f (po) + (po)(p - Po) (23) 
where df is the differential of f and is given as: 

It follows that in this neighbourhood the function f provides a 1 to 1 
relationship. Thus, if it is possible to establish that the Jacobian is strictly 
positive or strictly negative throughout the whole of P-space, the function f 
will be 1 to 1 everywhere. Once this condition has been established, it is 

10 necessary to consider how the error associated with image acquisition maps 
under f \ to give the corresponding error in parameter recovery. The error 
associated with each image ratio is calculated using equation (19). The 
combination of errors maps out a hypervolume in I -space, centred on the 
point io. This has been illustrated in Figure 2 for the case of a 2D P-space, 

15 where an ellipse is obtained (or a circle if the errors are equal). An 
ellipsoid is obtained in 3D space and a hyperellipse in higher dimensions. 
Although the following analysis will be based on a 2D P-space, the 
arguments are equally valid in higher dimensions. 



20 The ellipse in I -space represents all possible image ratio vectors which 
could correspond to a camera measurement i = io- It is assumed that the 
region of error lies within the neighbourhood of i = where the mapping 
function f can be approximated linearly. Thus, under the mapping r\ the 
ellipse in I -space maps directly to another ellipse in P-space. This new 

25 ellipse is centred on the point p = po and represents all possible parameter 
vectors which could be recovered from the vector of image ratios i = io- 
The error associated with parameter recovery is obtained by considering the 
worst case scenario: that is the point within the ellipse in P-space which is 
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at the maximum distance from the point p = po. This maximum distance 
must be calculated separately for each component, pt, of the parameter 
vector to obtain the error associated with recovery of each individual 
component. To calculate these errors it is necessary to consider how the 
5 ellipse is transformed under the mapping r\ which is linear provided the 
Jacobian is non-zero. 



Under a linear mapping the ellipse will be translated, scaled and rotated. 
The translation associated v^ith the linear mapping defines the point p = po 
10 which is mapped to from the point i = io. The two other transformations, 
scaling and rotation, are best understood by considering how a vector 
di = df in I -space, maps under to give a corresponding vector dp in P- 
space. The vector dp can be calculated from the inverse form of equation 
(24) which, in matrix form, is given as: 

dp = J'^di (25) 

15 

where J denotes the Jacobian matrix. Note that J'^ exists only if the 
Jacobian in non-zero. This must be the case if the 1 to 1 condition is to be 
satisfied. 



20 



25 



The vectors A and B correspond to the major and minor axes of the ellipse 
in I -space and are given as: 



B = 




(26) 



Under the mapping these vectors map to the vectors A' and B' which 
correspond to the major and minor axes of the ellipse in parameter space. 
Solving equation (25) for each of these vectors gives: 



A' = 




(27) 
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where Apf and are the components of the vector A' in the direction of 
pi and p2 respectively. Similarly, Ap^ and Ap^ are the components of the 
vector B' in the direction of pi and p2 respectively. To calculate the error in 
each component of the parameter vector it is necessary to consider the 

5 worst case scenario. It can be seen from Figure 2 that this corresponds to 
taking the maximum of Ap/^ and Ap^ as the error in pi and taking the 
maximum Apj ^^^d Ap^ as the error in p2. This error can be specified by a 
vector Ap and can be calculated for any given point in parameter space. 
With this measure of the accuracy of parameter recovery across the whole 

10 of parameter space, it is possible to define an optimisation criterion. This 
could simply be based on a sum of the errors at every point in P-space or 
could be chosen to favour accuracy of recovery of a subset of the original 
parameters. Once this optimisation criterion has been defined, a search can 
be used to find the optimum mapping function f. It is important to note 

15 that, although the above discussion is based on a 2D parameter space, the 
methodology is equally applicable to any JS'-dimensional parameter space. 

An algorithm for the implementation of the proposed methodology is given 
as follows: 

20 1. Establish a suitable search space from a parameterisation of 



mappings b and c. 



2, For a given mapping function f calculate the vector of image 
ratios for each point within a discretised parameter space. 



25 



3. For each point, check that the Jacobian is either strictly positive or 
strictly negative across the whole of parameter space. If this 
condition is held then compute the inverse of the Jacobian matrix. If 
not then return to step 1 and define a new mapping function f . 



4. Using equation (19), calculate the error associated with image 
acquisition for each vector of image ratios. 



30 



5. From the image ratio vector error calculate the maximum possible 
error in each component, pjk, of the parameter vector across the whole 
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of parameter space. 

6. Use the vector of parameter errors at each point within parameter 
space to measure the accuracy of parameter recovery. 

7. Repeat steps 2-6 with some optimisation technique which enables 
5 an optimum mapping function f to be determined. 

It is fairly straightforward to extend this methodology to the case in which 
N > K : that is, where there are more image ratios than parameter values. 
Initially every possible K-dimensional subspace of image ratios will need to 
be defined from the original N -dimensional space of image ratios. It will 

10 then be necessary to go through the above procedure for each potential 
subspace and obtain the vector of parameter errors at each point within 
parameter space. To achieve the maximum possible accuracy the best Ap 
must be selected at every location within parameter space. Thus every 
point in P-space will be linked to a specific image ratio combination. It 

15 will then be necessary to link every region of the original N -dimensional 
space of image ratios to the particular subspace of image ratios which 
should be used for parameter recovery. It is important to note that it is 
necessary to recover the whole parameter vector at each point i^, within a 
particular J^-dimensional subspace of image ratios. It is not possible to 

20 attempt to improve the accuracy of the system by recovering different 
components of the parameter vector from different i5r-dimensional 
subspaces of image ratios. This is mathematically invalid. 



The mapping function f is a composite function of three separate mappings. 

25 Although the first mapping a, from parameter space to wavelength space, is 
fixed for a given problem, mappings b and c can vary depending on the 
choice of optical filters and definition of image ratios. Thus, to define an 
appropriate search space it is necessary to parameterise mappings b and c. 
Mapping b, which represents image acquisition, is defined by the positive 

30 NxM matrix R^, given in equation (11). Typically this matrix will 
contain many elements and an appropriate parameterisation should be based 
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on typical filter response functions. For exaraple, the position of the central 
wavelength and a measure of width could be used to define a Gaussian 
shape. 

5 Parameterisation of the mapping function c will be fairly straightforward as 
there are only a limited number of ways of combining image values to 
produce independent image ratios. In some applications the form of this 
mapping may be fixed apriori. Thus, it will not increase the overall 
dimensionality of the search space. 

10 

An optimisation method should search the whole space of possible 
mappings using the optimisation criterion outlined in the previous section. 
One technique which is ideally suited to this type of search is a genetic 
algorithm, GA, {see T. Back and H. P. Schwefel, "An overview of 

15 evolutionary algorithms for parameter optimisation," Evolutionary 
Computation 1, 1-23 (1993)} as it is straightforward to define a fitness 
function which measures the accuracy of parameter recovery. Genetic 
algorithms have been shown to work well on a wide range of problems with 
objective functions that do not possess "nice" properties such as continuity, 

20 differentiability or satisfaction of the Lipschitz Condition {see L. Davis, 
The handbook of genetic algorithms (Van Nostrand Reingold, New York, 
1991), and D. Goldberg, Genetic algorithms in search, optimization and 
machine learning (Addi son-Wesley, London, 1989)}. 



25 The above techniques will now be further exemplified by considering their 
application to the analysis of a body component, in this case a normal skin 
composition. Firstly, the prediction of spectral reflectance is considered. 



In order to perform mapping a it is necessary to have either a mathematical 
30 model which can predict spectral reflectance for a given set of parameter 
values or some technique for measurement of the appropriate spectrum. For 
this application we use the mathematical model developed by Cotton and 
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Claridge {see S. D. Cotton and E Claridge, "Developing a predictive model 
of human skin colouring," Proc. of SPIE Med. Imag. 2708, 814-825 
(1996)}. With this model it is possible to predict the spectral reflectance 
for a given set of parameters. An outline of the model is now given. 

5 

Skin can be considered to be the four-layer structure depicted in Figure 3. 
A negligible amount of light is reflected from the surface of the skin, thus 
the surface term in equation (2) can be neglected. Although not absorbing 
any radiation, the stratum comeum scatters the incoming light in all 

10 directions. Light which penetrates this layer can thus be considered diffuse. 
In the epidermis light is absorbed by the pigment melanin. The absorption 
at each wavelength can be calculated using the Lambert-Beer law and will 
depend on the product of the melanin extinction coefficient and the pigment 
concentration. After passing through the epidermis the light is both 

15 scattered and absorbed by the papillary dermis. The absorption results from 
the presence of blood and scattering from the underlying collagen structure. 
The simple Kubelka- Munk light theory {P. Kubelka and F Munk, "Ein 
Beitrag zur Optik der Farbanstriche", "Z. Tech. Opt" 11, 593-611 (1931)} 
can be used to model the interaction of light with the papillary dermis as 

20 the necessary condition of diffuse incident illumination is satisfied. Any 
light which passes through the papillary dermis into the recticular dermis 
can be neglected as no significant backscattering occurs in this layer. 
Using this two-layer light transport model it is possible to obtain the 
remitted spectra for given concentrations of melanin and blood. A more 

25 detailed description of this model can be found in S. D. Cotton and E 
Claridge, "Developing a predictive model of human skin colouring", Proc, 
of SPDE Med. Imag. 2708, 814-825 (1996). 

For a given papillary dermal thickness, changes in melanin and blood 
30 characterise all histological variation and thus define a 2-D parameter space 
for healthy skin. To carry out the optimisation procedure described above it 
is necessary to discretise parameter space. This is done at equal intervals to 
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define 10x10 points, each of which corresponds to a spectrum generated by 
the mathematical model. For simplicity, concentration values will be 
denoted by a number between 1 and 10. Figures 4a and 4b show how the 
remitted spectrum changes as melanin and blood are varied respectively. 
5 With a change in melanin concentration, the intensity of the whole spectrum 
is seen to decrease, with a more pronounced change in the blue region. As 
the blood concentration is decreased the most significant reduction in 
intensity is observed in the green region, the resulting shape reflecting the 
two absorption maxima of oxyhaemoglobin, a blood bom pigment. 



To define a suitable search space it is necessary to parameterise the 



interference filter. This is modelled as a square profile with Gaussian 
decay at each side. Two parameters are required to specify this shape: the 
15 central wavelength and a full width half maximum (FWHM). Optimisation 
is carried out for three such filters, defining a 6-D search space. With three 
filters giving three image values, ir, I'a and 13, the only possible definition of 
image ratios, if we assume 1*1 = 13 is equivalent to 13 = ii, is given as: 

= — ^2 = V-^ 

20 In this instance the mapping c does not increase the dimensionality of the 
search space. 

The optimisation procedure was implemented following the algorithm given 
above. Initially the vector of image ratios was calculated for every point 
25 within the discretised parameter space. This was done using the 
mathematical model to perform mapping a, the parameterised form of 
matrix to perform mapping b and the equations (28) to perform 

mapping c. The derivative of each image ratio, with respect to each 
parameter, was obtained at each point within discretised parameter space 
30 using three-point finite difference approximations. The Jacobian matrix 



10 



mappings b and c. A parameterised form of b is chosen to define a typical 
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was then constructed at every point within paranaeter space, and providing 
its determinant was non-zero everywhere, the inverse calculated. If this 
condition was violated then a new mapping f was defined. The errors 
associated with image acquisition were then calculated using equation (19). 

5 The absolute value of the error in each image value will vary depending on 
the camera gain setting. Although this constant will not affect the mapping 
f, it must be estimated in order to calculate the effective camera error. For 
this application it was taken to be 0.78% of the maximum value of all the 
image values across parameter space. This corresponds to a camera which 

10 has been set to give a maximum reading for the largest spectral reflectance 
and a camera error of two grey scale levels in an 8-bit representation. 

Using the procedure outlined above the error associated with parameter 
recovery in both melanin and blood was obtained for each point within the 

15 discretised parameter space. In order to find an optimum f, it is necessary to 
minimise the errors in recovery of both melanin and blood across the whole 
of parameter space. Thus the fitness function for the GA was taken to be 
the sum of the errors in both melanin and blood. This procedure was 
implemented in matlab™ using a standard GA to search the space of 

20 available mappings. 

The boundaries of the search space were chosen such that the central 
wavelength was constrained to lie in the visible region (400nm-700nm) and 
such that the widths of the filters were allowed to vary from a FWHM of 25 
25 to 200nm. Although it is now possible to engineer almost any shape of 
interference filter, this corresponds to an economically viable range of such 
filters. 

Although it was originally assumed that an image ratio defined as ii I h 
30 would be equivalent to 13 / z'l, the results of the GA search showed that this 
was not the case. The search was intitialised for a random seed and, 
although the same central wavelengths were always obtained, different 
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filters were selected corresponding to 13 defined in equation (28). Further 
investigation showed that these local maxima in the search space 
corresponded to differing distributions of errors both, across parameter 
space and between the two parameters. This is because the fitness function, 

5 or measure of accuracy, was defined as the sum the errors across parameter 
space for both melanin and blood. Thus, a loss of accuracy in one 
parameter could be compensated for with an increase in the other. It may 
be that, with a more exact specification of the error distribution in the 
fitness function, it would be possible to obtain the same results for every 

10 GA search. 

Figure 5 shows a filter combination which gave a similar error in the 
recovery of both melanin and blood. The image ratios were calculated by 
dividing the filter centred at X = 473nm and X = 560nm by the response of 

15 the filter centred at X = 700nm. To understand why these specific filters 
were selected it is necessary to analyse the spectral curves shown in Figure 
4. The filters centred at A. = 473nm and X = 560nnmm correspond to 
spectral locations where there is a large change in intensity with the 
parameters melanin and blood respectively. A third filter was then required 

20 in a region of the spectrum in which the remitted light which was either 
significantly less or significantly more than that of the other two filters. 
The filter centred at X = 700nm was chosen as it always gave the largest 
response at any point within parameter space. This ensured that the 
derivatives of each image ratio decreased monotonically across the whole of 

25 parameter space. The Jacobian, calculated from these derivatives, was 
strictly positive across the whole of parameter space. It is interesting to 
note that some alternative filter combinations gave Jacobians which were 
strictly negative across parameter space, corresponding to alternative local 
maxima in the search space. If two filters are chosen, to define an image 

30 ratio, which vary similarly across parameter space, there will be minimal 
change in that image ratio and thus it will be of limited value for parameter 
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measurement. 

It has been demonstrated that, using an objectively defined set of optical 
filters, it is possible to recover scene parameters from image data in a way 

5 which is insensitive to geometry and incident illumination. In the example 
problem, discussed above, the error associated with this parameter recovery 
was found to be relatively small. The invariance of this mapping means 
that the technique will be particularly applicable to medical image 
applications where there is significant curvature of the surface of the 

10 imaged tissue, such as near a joint. It also means that the method can be 
used for whole body imaging. It will also be unnecessary to calibrate the 
camera to determine the intensity of the incident light. This could help to 
significantly increase the speed of image acquisition and later processing. 

15 The methodology set out here has been developed for a measurement task, 
where the scene parameters are known to vary continuously. The technique 
can be also be applied to problems of recognition, where it is necessary to 
differentiate discrete objects based on some measure of their spectral 
reflectance. This approach has been discussed in the article G. Healey, 

20 "Using colour for geometry-insensitive segmentation," J. Opt. Soc. Am. A 
6, 920-937 (1989) who used the idea of normalised colour to identify 
different regions of normalised colour space corresponding to different 
metal and dielectric materials. This enabled geometry-insensitive 
segmentation of an imaged object comprised of a number of different 

25 materials. 

It will be appreciated that in order to implement the proposed methodology, 
a look-up table should be established between all possible image ratios and 
scene parameters. Although this may be time consuming, it is only 
30 necessary to carry out this procedure once. Once established, this look-up 
table will ensure no significant processing after image acquisition, making 
this technique particularly suitable to real-time applications. 
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10 



15 



20 



Figure 6 is a flow diagram showing the key steps in the method described 
above. 

The method and apparatus for analysing at least one parameter of a body 
component, in this case animal tissue in the specific form of facial skin is 
illustrated in figure 8. A light source 100 provides illumination to the tissue 
and remitted light is received at photoreceptor 200 which in this case is a 
digital camera. Two cross polarised linear polarising filters 300 are used to 
eliminate the effects of surface reflection from the skin. One filter 300 is 
placed between the light source 100 and the skin and the other filter 300 is 
placed between the skin and the digital camera 200. 

In this case the digital camera is provided with Red, Green and Blue filters 
so that light in those wavebands is received by the camera. These 
wavebands are used to provide image ratios of which the concentration of 
melanin and the concentration of blood are one to one function. 

The procedure outlined in figure 7 was applied to image data in the 
following way. 

1. Two parameters: the concentration of melanin and blood were 
identified as sufficient to describe all histological variation of 
healthy tissue. 

2. A Kubelka-Munk model of light transport was used to predict the 
remitted spectrum of tissue for any given combination of melanin 
and blood concentration. 

3. The spectral responses of each of the RGB channels of the colour 
camera were established and image ratios defined as 



a. Ratio 1 = Green / Re <i 



Ratio 2 = Blue/Red 
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4. 



The mapping from the 2-D space of parameter variation to the 2-D 
space of image ratios was checked to ensure that is was 1-1 across 
the whole range of appropriate parameter variation. 



5 



5, 



A piecewise continuous approximation was constructed to define a 
function relating image ratios to histological parameters. 



6. 



Images were acquired using a system of crossed polarising filters, as 
described above. The experimental set up has been illustrated in 
figure 8. 



10 



7, 



The function described in step 5 was then used to process the image 
data. 



8. Parametric maps were then produced of melanin and blood across the 
imaged tissue. 

In one experiment this method was applied to an image obtained using a 
15 JAI CV-M7CL+ camera imaging facial skin. Parametric maps, showing a 
grey-scale representation of then quantitative measurements of melanin and 
blood derived using this technique, are shown in figure 9b and 9c. 

It should be noted that in 9b illustrating the concentration of haemoglobin 
20 concentration across the image, spot S is identified but mole M is not 
identified. However in 9c illustrating the concentration of melanin across 
the image, spot S is not identified while mole M is identified. This 
illustrates simply how useful a tool this can be for a clinician. 

25 A second specific embodiment involves the analysis of images of the human 
gastrointestinal track obtained using an endoscope. The endoscope system 
can take two alternative embodiments. In one case the endoscope is 
equipped with a conventional colour camera and white light source 
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equipped with cross polarizing filters^ In a second case the endoscope is 
equipped with a monochrome camera and a light source equipped with cross 
polarizing filters^ with the light source that changes colour sequentially 
between red, green and blue, and these changes are synchronised with the 
camera to produce a sequence of red, green and blue images. 



The procedure outlined in figure 7 is applied to this problem, using data 
from an endoscope equipped with a conventional colour camera, as follows: 

10 1. Two parameters blood concentration and tissue thickness are 

identified as sufficient to describe all histological variation. 
2. A Monte Carlo model of light transport is used to predict the 
remitted spectrum of the given tissue for any possible combination of 
blood concentration and tissue thickness. 
15 3. For an endoscope and camera system, the spectral responses of each 

of the RGB channels is established and image ratios defined as 
a. Ratio 1 = Green /Red Ratio 2 = Blue /Red 

4. The mapping from the 2-D space of parameter variation to the 2-D 
space of image ratios is checked to ensure that it is 1-1 across the 

20 whole range of appropriate parameter variation. 

5. A piece wise continuous approximation is constructed to define a 
function relating image ratios to histological parameters. 

6. Images are acquired the endoscope with a system of crossed 
polarising filters. 

25 7. The function described in step 5 is then used to process the image 

data. 

8. Parametric maps are then produced to display variation in blood and 
tissue thickness across the given image. 



' one filter being placed between the source of illumination and the component, and the 
other filter placed between the component and the photoreceptor or photoreceptors with 
the filters being set at 90 degrees to one another. 
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The procedure can be modified to analyse additional histological 
parameters with the addition of additional wavebands as described in the 
equations shown above. These additional wavebands may be obtained by a 
5 monochrome camera and light source with cross polarising filters taking a 
series of images of the subject illuminated by a sequence of coloured lights 
of known spectral characteristics. The spectral characteristics of one or 
more of colours may lie outside the visible spectrum. 
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Claims 



1. A method of analysing at least one parameter of a body component, 
comprising the steps of illuminating the component or body with light of at 

5 least a first and second waveband, receiving light of at least said first and 
second wavebands remitted by the component at a photoreceptor or 
photoreceptors, and analysing the light received at the photoreceptor(s) to 
provide a ratio between the amount of light remitted of the first waveband 
and the amount of light remitted of the second waveband, and from this 

10 calculating the component parameter. 

2. A method according to claim 1, in which the wavebands are 
predetermined and calculated by use of a mathematical model of the body 
component and its characterising parameters. 

15 

3. A method according to claim 1, in which the wavebands are 
predetermined and derived through use of a biological model of the body 
component. 

20 4 A method of analysing at least one parameter of a body component, 
comprising the steps of illuminating the body or component with light of at 
least a first and second waveband, receiving light of at least said first and 
second wavebands remitted by the component at the photoreceptor(s), but 
eliminating light reflected by the component or body and analysing the light 

25 received at the photoreceptor(s) to provide a ratio between the amount of 
light of the first waveband and the amount of light of the second waveband, 
and from this calculating the component parameter. 

5 A method according to any of the preceding claims where the 
30 wavebands are chosen such that the component parameter is a one to one 
function of the ratio between the amount of light remitted by the body 
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component of the first waveband and the amount of light remitted by the 
component of the second waveband. 

6 A method according to any of the preceding claims, in which the waveband 
5 ratios are compared with a mathematically generated model of waveband ratios 

corresponding to a range of component parameters. 

7 A method according to any of the preceding claims, in which the waveband 
ratios are compared with an experimentally measmred set of waveband ratios 

10 corresponding to a range of component parameters. 

8 A method according to claims 6 or 7 where the comparison results in a 
measure or measures relating to the component parameter or parameters. 

15 9 A method according to any of the preceding claims where a function is 
derived relating the computed ratios and the component parameter or parameters. 

10 A method according to any of the preceding claims, in which the 
light reflected by the component is eliminated by the use of a pair of cross 
20 polarised linear polarizing filters, one filter being placed between the 
source of illumination and the component, and the. other filter placed 
between the component and the photoreceptor or photoreceptors. 

11. A method according to any of the preceding claims in which the light 
25 illuminating the body component is a light of a plurality of wavelengths 
which includes at least the wavebands. 

12 A method according to claim 11 in which the illuminating light is 
ambient light. 



30 



13 A method according to claim 11 in which the illuminating light is 
sunlight. 
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14. A method according to claims 11 - 13, in which at least one filter is 
placed sequentially between the source of illumination and the component 
or between the component and the photoreceptor or photoreceptors. 

5 

15. A method according to any one of the preceding claims, in which the 
body component is human or animal tissue. 

16. A method according to claim 15, in which the tissue is one of skin, 
10 the lining of the gut, colon, oesophagus, cervix, eye or any other epithelial 

tissue. 

17. A method according to any one of the preceding claims, for 
analysing a plurality of body component parameters where the component is 

15 illuminated with light of each of a collection of wavebands. The light 
received by the photoreceptor or photoreceptors includes this collection of 
wavebands. The light is analysed at the photoreceptor(s) to provide a 
collection of ratios between the amount of light of each waveband with 
some or all of the other wavebands and from this calculating the component 

20 parameter(s). 

18. A method according to any one of the preceding claims, for 
analysing a plurality of body component parameters in which for each 
component parameter there exists a pair of predetermined wavebands such 

25 that the component parameter is a one to one function of the ratio between 
the amount of light remitted by the component of the first predetermined 
waveband of the pair and the amount of light renodtted by the component of 
the second predetermined waveband of the pair, and the component is 
illuminated with light of each pair of predetermined wavebands, the light 

30 received by the photoreceptor or photoreceptors is of each pair of 
predetermined wavebands remitted by the component at the 
photoreceptor(s), and analysing the light received at the photoreceptor(s) to 



wo 2004/010862 ^^CT/GB2003/003367 

44 

provide for each component parameter a ratio between the amount of light 
of the first waveband and the amount of light of the second waveband, and 
from this calculating each component parameter. 



19 A method according to claim 17 or 18, in which the waveband ratios are 
compared with a mathematically generated model of waveband ratios 
corresponding to a range of component parameters to ascertain the component 
values. 

20 A method according to claim 17 or 18, in which the waveband ratios are 
compared with an experimentally measured set of waveband ratios corresponding 
to a range of component pairameters to ascertain the component values. 

21 A mefliod according to claims 19 or 20 where the comparison results in a 
measure or measures relating to the component parameter or parameters. 

22 A method according to claims 19 or 20 where a function is derived 
relating the computed ratios and the component parameter or parameters. 

23. A method according to any of the preceding claims in which the 
body component is skin and the parameters are the concentration of melanin 
and the concentration of blood. 

24. A method according to claim 23, in which the predetermined 
wavebands are the Red, Green and Blue colour bands, with the three 
wavebands providing the two ratios which are a one to one function with 
the parameters. 

25. A method according to any one of the preceding claims, in which the 
predetermined wavebands have been calculated by the steps of: 
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1) defining a set of potential wavebands 

2) defining one or more image ratios, the or each image ratio for a 
region being obtained by dividing the amount of light remitted by the 
component of a given waveband for that region, the "image value" 

5 for that filter, by another image value for that same region; 

3) for a component parameter to be analysed and for said defined set 
of potential wavebands and for said given image ratios, obtaining a function 
mapping points in parameter space to points in image ratio space; 

4) determining whether the mapping function provides a 1:1 
10 correspondence between points in parameter space and points in image ratio 

space; and 

5) if the mapping function does not provide a 1:1 correspondence, 
rejecting said potential wavebands, repeating steps 1) to 4) and, if the 
mapping function does provide a 1:1 correspondence accepting the potential 

15 wavebands as a candidate set of predetermined wavebands. 

26. Apparatus for analyzing at least one parameter of component, 
comprising a light source for illuminating the component with light of at 
least a first and second predetermined waveband, a photoreceptor or 

20 photoreceptors for receiving light of at least said first and second 
predetermined wavebands remitted by the component reflected by the 
surface at a photoreceptor or photoreceptors; surface reflection elimination 
means for eliminating light reflected by the surface of the component and 
means for analyzing the light received at the photoreceptor(s) to provide a 

25 ratio between the light of the first waveband and the light of the second 
waveband , and from this calculating the component parameter. 



27 A method according to claim 26 where the predetermined wavebands 
are chosen such that the component parameter is a one to one function of 
30 the ratio between the amount of light remitted by the component of the first 
predetermined waveband and the amount of light remitted by the component 
of the second predetermined waveband. 
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28. Apparatus according to claim 26, in which the photoreceptor 
comprises a digital camera. 

29. Apparatus according to claim 28, in which the digital camera 
includes a plurality of filters, one for each predetermined waveband. 

30. Apparatus according to claim 26 in which the light source is ambient 
light. 

31. Apparatus according to claim 27, 28, 29 or 30, in which the distance 
between the photoreceptor(s) and the component is between 0.5cm and 
10m. 

32. Apparatus according to any one of claims 33 to 34, in which the 
distance between the light source and the component is between 0.5cm and 
10m. 

33. A method for deriving a pair of predetermined wavebands suitable 
for use in analysing a given parameter of a body component, the method 
comprising the steps of: 

1) defining a set of potential wavebands 

2) defining one or more image ratios, the or each image ratio for a 
region being obtained by dividing the amount of light remitted by the 
component of a given waveband for that region, the "image value" for that 
filter, by another image value for that same region; 

3) for the parameter of the component to be analysed and for said 
defined set of potential wavebands and for said given image ratios, 
obtaining a function mapping points in parameter space to points in image 
ratio space; 
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4) determining whether the mapping function provides a 1:1 
correspondence between points in parameter space and points in image ratio 
space; and 

5) if the mapping function does not provide a 1:1 correspondence, 
5 rejecting said potential wavebands, repeating steps 1) to 4) and, if the 

mapping function does provide a 1 : 1 correspondence accepting the potential 
wavebands as a candidate set of predetermined wavebands. 

34. A method according to claim 33 and comprising, for a set of 
10 potential wavebands accepted as a candidate set of wavebands, determining 
the accuracy of parameter recovery obtained using said mapping function 
and determining whether or not the accuracy is sufficient or matches some 
other criterion. 

15 35. A method according to claim 34, wherein if the accuracy is sufficient 
or matches said criterion, the candidate wavebands are adopted and if the 
accuracy is not sufficient or does not match said criterion, steps 1) to 5) are 
repeated for a different set of wavebands. 

20 36. A method according to claim 34 or 35, wherein the accuracy is 
determined by: 

a) calculating the error associated with image acquisition for each 
vector of each image ratio; 

b) from the image ratio vector error, calculating the maximum 
25 possible error in each component of the parameter vector across the 

whole of parameter space; and 

c) using the vector of parameter errors at each point within 
parameter space to measure the accuracy of parameter recovery. 



30 37. A method according to claim 33 or 34 and comprising repeating 
steps 1) to 5) for a multiplicity of sets of potential wavebands to identify a 
plurality of candidate waveband sets, determining for each candidate set an 
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error value representing the accuracy of parameter recovery obtained using 
the corresponding mapping function, and using said candidate set as a basis 
for determining a preferred set of wavebands. 



5 38 A method according to claim 37 and comprising using a genetic 
algorithm to determine a preferred set of wavebands using said candidate 
set. 



39. A method according to any one of claims 33 to 35 and comprising 
10 using a gradient descent algorithm to select an optimal set of wavebands, 

the starting point for the algorithm being a first candidate set of wavebands 
identified in step 4). 

40. A method according to any one of claims 33 to 39, wherein for each 
15 waveband the method of the present invention is used to determine the 

center wavelength of the waveband. 



41. A method according to any one of claims 33 to 40, wherein the 
method is used to determine the full width half maximum (FWHM) of the 
20 waveband. 



42. A method according to any one of claims 33 to 41, wherein step 3) 
comprises constructing a Jacobian matrix for the mapping function with 
respect to said parameter(s), and obtaining the determinant of that matrix. 

25 

43. A method of determining a property or properties of each of a set of 
filters, which filters are used to select specific wavelength ranges in a 
system which relies upon a spectral analysis of remitted, emitted, and/or 
transmitted light to quantify a parameter or parameters of an object or 

30 material, the method comprising the steps of: 

1) defining a set of potential filter properties 
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2) defining one or more image quotients, the or each image quotient 
for a region being obtained by dividing the quantified output of a given 
filter for that region, the "image value" for that filter, by another image 
value for that same region; 



potential filter properties and for said given image quotients, obtaining a 
function mapping points in parameter space to points in image quotient 
space; 

4) determining whether the mapping function provides a 1:1 
10 correspondence between points in parameter space and points in image 

quotient space; and 

5) if the mapping function does not provide a 1:1 correspondence, 
rejecting said potential filter properties, repeating steps 1) to 4) and, if the 
mapping function does provide a 1:1 correspondence accepting the potential 

15 filter properties as a candidate set of filter properties. 

44. Apparatus for analysing an object or material having means for 
conducting a spectral analysis of remitted, emitted, and/or transmitted light 
to quantify a parameter or parameters of an object or material, the apparatus 
20 comprising a plurality of filters for splitting said light into respective 
components, the filters having properties obtained by using the method of 
the claim 43. 



5 



3) for an object or material to be analysed and for said defined set of 
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Figure 5 



wo 2004/010862 




CT/GB2003/003367 



Define discrete representation of spaces P, A, I and t 
and define a mapping function a: P — > A 



Define filter properties (a mapping function b: A I) and 
define image quotients (a mapping function c: I — > I) 



^ 

For the given definition of filter properties and the given 
definition of image quotients obtain a mapping function 

f:P-^ t 



Using the given mapping function f calculate the vector of 
image quotients for each point (parameter vector) p within 
the discretised parameter space P 



For each point p in P compute the Jacobian d^dp 



^ Determinant of the ^ 
Jacobian strictly positive or strictly 
^negative for each point in P 



Amend filter properties 
and/or image quotients 



YES 



For each vector of image quotients calculate the error 
associated with image acquisition 



From the image quotient error calculate the maximum error in 
each component pk of the parameter vector p across the whole 
of parameter space P 

i 



Measure the accuracy of parameter recovery using the vector 
of parameter errors at each point p in parameter space P 



Accuracy of 
parameter recovery optimal and/or 



Amend filter properties 


^ NO ^^^^ 


and/or image quotients 





YES 



Figure 6 
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For the tissue to be analysed identify all histological 
parameters which could cause a variation in the 
spectrum of remitted light. 



Using an experimental and/or theoretical technique have 
some means of predicting a remitted spectrum for a 
given combination of histological parameters. 

Establish the spectral response of each of the 
channels/filters used in image acquisition system 





Obtain an image of the scene. From a consideration of 
image values recorded by the camera, define some 
suitable image ratios 






Have some means of calculating the Jacobian matrix for 
any point in parameter space 







Define alternative image ratios 

and/or choose an image 
acquisition system with different 
spectral responses 



Determinant non-zero 
across the whole of 
parameter space 



NO 



YES 



Construct a function relating every possible image ratio 
to every possible parameter combination. 



Construct a function relating every possible image ratio 
to every possible parameter combination. 



Use this function to process image data to give 
quantitative information on histological parameters 



Figure 7 
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Figure 8 
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